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Abstract 

This paper presents a lattice approach to model the influence of cracking on inviscid flow in concrete. 
A mechanical lattice model based on a damage-plasticity constitutive model was combined with a new 
dual lattice of conduit elements for flow analysis. The diffusivity of the conduit elements depends on the 
crack-opening obtained from the mechanical lattice. The coupled lattice model was applied to several 
benchmark tests for aligned and random lattices. The results for mechanical loading and flow analysis 
obtained with the new approach were shown to be independent of the size of lattice elements used. 
Keywords: Concrete, cracks, lattice, transport, flow, Voronoi tesselation 



1 Introduction 

Cracking increases the permeability and diffusivity of concrete, which may accelerate the deterioration 
of concrete structures 31, 11, [29|, [l6j]. Cracks can act, for instance, as pathways for water containing 



chlorides, which promote corrosion of reinforcement and subsequently may lead to further cracking and 
spalling of the concrete cover. For cracking induced by mechanical loading, macroscopic crack widths 
in concrete are a result of complex fracture processes on the meso-scale. For instance, crack bridging 
and crack tortuosity depend on the material composition and may influence permeability. Other types 
of loading, such as drying shrinkage subjected to aggregate restraint, lead to distributed micro-cracking, 
which depends on aggregate size and volume fraction and cannot be uniquely related to a macroscopic 
crack width 32, Numerical analysis of the fracture process of concrete on the meso-scale within 



the framework of the finite element method may increase the understanding of the influence of material 
composition on micro-cracking and transport properties. However, this requires the development of robust 
and accurate analysis methods for the interaction of fracture and mass transport. 

Modelling of flow in cracked concrete within the finite clement framework can be categorised according 
to nonlinear fracture mechanics approaches. In continuum approaches, fracture is represented as regular 
fields of localised inelastic strains of finite size by using higher order constitutive laws, such as integral- 
type nonlocal models 0, Q ■ Transport processes in the continuum can be spatially varied according to 
the values of accumulated history variables, such as inelastic strains. This approach is computationally 
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demanding since a fine discretisation is required to be able to describe the localised strain fields. Al- 
ternatively, the representation of the fracture process zone is simplified to a displacement discontinuity 
embedded into the continuum [H El 0] ■ For these hybrid approaches, localised transport along the dis- 



crete crack is coupled with continuum modelling of imbibition into the adjacent material [2a. l8l l28l. l26l] . 
Coupling of the flow in discrete cracks and the surrounding material can create difficulties in these ap- 
proaches [2(| . Finally, in discrete approaches the fracture process is described by the failure of structural 
elements, such as trusses and beams. Lattice models have shown to be capable of describing complex 
fracture patterns on the meso-scale of concrete 15, 27 , Hoj |. Furthermore, mass transport can be de- 



scribed by a lattice of conduit elements, which can be linked to the structural lattice to couple fracture 
and transport processes [9]. Some lattice approaches exhibit mesh dependency and are limited in de- 
scribing the continuum response accurately. A special type of lattice model for the mechanical response 
and mass transport was proposed recently, which provides mesh-independent and accurate descriptions 
of basic aspects of the continuum response @, @, 0| • In this approach, the cross-sections of structural and 
transport elements are determined from the Voronoi tessellation of random points placed in the strutural 
domain. This modelling approach was further developed to describe the interaction of transport along 
discrete cracks and the surrounding material by introducing an additional lattice of transport elements 



21|, |30(. This research direction is analogous to the hybrid approaches, with the challenges arising in the 



interaction of the two lattices describing transport along cracks and transport in the adjacent material. 

In the present study, a new lattice approach is developed to provide a mesh-independent description of 
the fracture process and mass transport in the cracked material. Transport is modelled by one lattice, 
which is dual to the lattice used for the mechanical response. The change of transport properties due 
to localised crack openings is smeared out over the width represented by the lattice element, ensuring a 
mesh-independent description. The first objective of this work is to demonstrate that a lattice dual to 
the mechanical network accurately represents mass transport in the uncracked continuum. The second 
objective is to investigate if flow in a cracked media, such as concrete, can be represented independently 
of the size of the element mesh. To the author's knowledge, this is the first lattice approach that models 
flow in fractured media independently of the mesh size. 



2 Combined model for mechanical loading and flow 

2.1 Mechanical model 

In this work the mechanical response of concrete is described by a lattice approach [3, [j| with a damage- 
plasticity constitutive model [12j . In the following sections the modelling approach is briefly reviewed. 

The lattice model is based on the Voronoi tesselation of the domain [22|. Lattice elements connect nodes, 
which are the nuclei of the neighboring Voronoi polygons (Fig. [T£i) . The nodes are placed randomly in the 
domain, constrained by a minimum allowable distance d m between nodes. The smaller d mi the smaller 
is the average distance between the nodes and the average size of the lattice elements Each node 
has three degrees of freedom (Fig. [T]d) , that is two translations and one rotation, which determine the 
displacement jump at midpoint C of the element cross-section in the local coordinate system as 

u c = Bu (1) 

where 

[-1 e 1 -e 1 

[ -1 -h/2 1 -h/2\ ( ' 

and u c = {mi, ui, 0i, «2, f2, 02} T - In Eq. (3J e is the eccentricity of the midpoint C with respect to the 
element axis and h is the length of the element (Fig. [TJd). If C is on the left hand side of the element, 
e = e c . Otherwise, e = — e c . 
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Figure 1 : Discretisation: (a) Elements (solid lines) and corresponding cross-sections (dashed lines) of the 
mechanical lattice obtained from the Delaunay triangulation and dual Voronoi tessellation, respectively, 
(b) Degrees of freedoms of the lattice elements for the mechanical and transport model in the local 
co-ordinate system. 



The cross-sections of the lattice elements are the faces of the Voronoi polygons. The element stiffness in 
the local coordinate system is 

K c = ^-B T D C B (3) 

h 

where I is the width of the cross-section (polygon facet) and D e represents the material properties as 
described below. The present model was proposed for two-dimensional plane-stress analysis; for the 
examples considered later, the out-of- plane thickness was assumed to be equal to 1. 

The constitutive model, which relates the strain vector e = u c /h to the nominal stress vector er, is based 
on a combination of plasticity formulated in the effective stress space and isotropic damage mechanics. 
The stress-strain law is 

<T = (l-w)D e (e-£ p -e T ) = (l-w)<r (4) 

where ui is the damage variable, D is the elastic stiffness, e p = (e pn ,£ps) T is the plastic strain, £t = 
(st,0) T is the eigenstrain, and a — (a D ,cr s ) T is the effective stress. The subscripts n and s denote the 
normal and shear direction in the local coordinate system (Fig. []J>) . The eigenstrain is defined here as a 
strain which is not due to mechanical loading, but to other physical or chemical processes. The elastic 
stiffness is 

Mo :*} < 5 > 

where E and 7 are model parameters controlling both the Young's modulus and Poisson's ratio of the 



material [14(. For instance, for plane stress considered here and a lattice of equilateral triangles, Poisson's 
ratio v is 

(6) 

3 + 7 

The plasticity part of the damage-plasticity model is based on the effective stress and consists of the 
yield surface, flow rule, evolution law for the hardening parameter, and loading-unloading conditions. 
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The yield surface is elliptic and its initial size and shape are determined by the tensile strength / t , the 
shear strength sf% and the compressive strength eft of the material under consideration. The evolution 
of the yield surface during hardening is controlled by the model parameter fi, which is defined as the 
ratio of permanent and total inelastic displacements. The scalar damage part is chosen so that linear 
stress inelastic displacement laws for pure tension and compression are obtained, which are characterised 
by the fracture energies Gft and G{ c . The equivalent crack opening is defined as w c = ||w c ||, where 

w c = h (e p + uj (e - e p )) (7) 

The crack opening vector w c is composed of a permanent and reversible part, defined as he p and 
uih (e — e p ), respectively [12 1. 



The constitutive behavior of the damage-plasticity model is illustrated by its stress-strain response for 
fluctuating normal strains for /i = 1 and /i = (Figure [5]) . The normal strain is increased to point A for 




Figure 2: Stress-strain response for fluctuating normal strains for \i = 1 (solid line) and /i = (dashed 
line) . 

the case of /i = 1 (and to point A for /j, = 0). A and A coincide, since the response for both parameter 
settings has been identical up to this point. The normal strain is then reduced to point B (B ) and again 
increased to point C (C ). For \i = 1 the combined approach reduces to a pure plasticity model. The 
unloading is elastic and under subsequent compressive loading the compressive strength is reached sooner 
than for /i = 0. On the other hand, a pure damage-mechanics response is obtained for \i — 0. The stress- 
strain curve is unloaded to the origin. This constitutive model in combination with the present lattice 
approach results in a mesh-independent description of the mechanical response. A detailed description 
of the components of the model is presented in [13] • 



2.2 Coupling of flow and mechanical models 

Here, flow is modelled by a new lattice approach, for which the spatial arrangement of the conduit 
elements is dual to the arrangement of structural elements. The conduit elements are placed along the 
facets of the Voronoi polygons (Fig. [T^,) . Equivalently, cross-sections of conduit elements are positioned 
on the mechanical lattice elements. The discrete form of the differential equation of the non-stationary 
flow problem for one conduit element is defined as 

BO 

a e 6 + C c — =f (8) 
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where a c and C c are the element conductivity and capacity matrix, respectively, t is the time, f are the 
external fluxes, and the degrees of freedom of the conduit elements are the flow potential 6 = (0i,62) T 
(Fig-El 3 ) El- The conductivity matrix is defined as 



ae= r(-i i 1 ) (9) 

where h is cross-section width, I is the length of the pipe element and a is the conductivity of the material. 
The capacity matrix C c is 

hi (2 1 



C *=12 1 2 



Mechanical loading is assumed to influence the diffusivity of the conduit elements as 

a = a + a c (h) (11) 

where ao is the initial diffusivity of the undamaged material and a c is the change of diffusivity due to 
mechanical loading. The part a c differs strongly depending on the problem modelled. For moisture 
transport, for instance, it could be related to the cubic law [31]. In the present study a simple linear law 
of the form 

a c = a - — (12) 
fi£fk 

was chosen, where w c is the equivalent crack opening from Eq. ([7|) and £fk is a parameter which controls 
the slope of the change of diffusivity. The equivalent crack opening is determined in Eq. ([7]) from the 
structural lattice element which crosses the conduit element. Since the structural and transport lattices 
are dual, the cross-section width h of the conduit element is equal to the length of the structural lattice 
element. An important aspect of the proposed model is the mesh independence, which is achieved by 
introducing the cross-section h in the expression of the diffusivity in Eq. (|12p . so that the change of 
conductivity due to mechanical loading in Eq. ([5]) is independent of the width of the conduit elements 
used. The parameters ao and £fk are material parameters. 



3 Model results 



The modelling approach described above, which was implemented in the object oriented finite element 
code OOFEM 23, 24 1, was applied to three benchmark problems. In the first example, a stationary flow 
field is represented on a random lattice. The second example involves the coupling of mechanical loading 
and flow for a lattice, which is aligned to a potential crack path. Finally, mechanical loading and flow are 
described for radial cracking about an inclusion within a random lattice. The numerical results of the first 
example are compared to the analytical solution. For the other two examples, possible mesh-dependence 
of the numerical results is investigated. 



3.1 Stationary flow within homogeneous media 

A graded lattice with conduit elements on the facets of Voronoi polygons is shown in Fig. 3. The cross- 
sectional areas of the conduit elements were chosen in two ways. In the first approach, cross-sectional 
areas were determined by the dual Delaunay triangulation as described in Section 12.21 In the second 
approach, a constant cross-sectional was used, which is the average of the areas obtained from the first 
approach. The nodes on the left and right hand sides of the model domain of length L = 0.1 m were 
subjected to constant potentials of 6 = and 6 — 1, respectively. For the other two edges, the boundary 
flux was assumed to be zero (q n — 0). The diffusivity was chosen as a = ao = 1 m 2 /s. The exact solution 
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Figure 3: Discretisation of the domain by means of a graded network. Lattice elements are located on 
the facets of the Voronoi polygons. 

for this problem is 9 = 10a;. The flow along the x-direction for y = L/2 for the two approaches is shown 
in Fig. UJ Additionally, the accuracy of the modelling approach was assessed by comparing the Li error 
norm to the exact solution. The error norm is 



and x is the position of node /, and 9 and 9 h are the exact and numerical values of the potential, 
respectively. Furthermore, N is the number of nodes in the specimen. The error for a constant cross- 
sectional area is e r = 0.052, whereas the error for cross-sectional areas obtained from the Delaunay 
triangulation is e r = 5.795 x 10 -10 . 

Consequently, the lattice with transport elements placed on the edges of Voronoi tesselation and cross- 
sections obtained from the dual Delaunay triangulation results in an accurate description of the stationary 
flow field. The results of the present study complement results obtained from a dual lattice approach, in 
which the conduit elements are placed on the edges of the Delaunay triangulation and the cross-sections 
are the facets of the Voronoi tessellation [Jl. 



3.2 Nonstationary flow along a planar crack 

The second example involves the coupling of fracture and nonstationary flow for an aligned lattice. The 
analysis is divided into two steps. In the first step, a square specimen of length L = fl.lm was subjected to 
an eccentrically applied tensile force F. The eccentricity was chosen as L/4 with respect to the centerline 
of the model domain (Fig. The parameters for the mechanical model were set to E = 40 GPa, 




(13) 



where 




(14) 



G 




x-coordinate [m] 

Figure 4: Comparison of the potential 9 along the cc-direction (y = L/2) for cross-sections obtained with 
the Delaunay triangulation and a constant average cross-section. 

7 = 0.33, f t = 4 MPa, G ft = 100 N/m, q = 2, c = 10, G fc = 50000 N/m. The material parameters 
result in a mechanical response, which is typical for concrete. Three mechanical lattices with minimum 
nodal distances of d m \ n = 0.008, 0.004 and 0.002 mm were chosen, with the medium density mesh shown 
in Fig. [5^. The mesh was aligned in the middle of the specimen, so that the lattice elements were 
perpendicular to the crack path. The damage-plasticity model was used only for elements crossing the 
predefined crack path. All other elements were assumed to be elastic. This allows one to evaluate the 
flow potential along the crack. The result in the form of the load F versus the displacement d at the 
loading point is shown in Fig. for three lattices. The load-displacement curves are almost completely 
independent of the size of the lattice elements. At a displacement of d = 0.1 m (marked in Fig. E]by an 
open circle), the crack almost reaches the right side of the specimen. For this stage, nonstationary flow 
analyses were performed for the three lattices dual to the ones used for the mechanical analyses. The 
lattice with the medium element size is shown in Fig.jSJx The nodes on the left hand side of the specimen 
were subjected to 9 — 1 at all times, whereas the other nodes have an initial potential of 9 = at t = 0. 
The diffusivity a was chosen according to Eq. (fT2"|) with ao = 1 m 2 /s and £fk = 0.0025. The results are 
presented by means of the potential 9 in x-dircction along the crack path (y — 0.05 m) in Fig. [7] and 
the potential 9 in y-direction perpendicular to the crack at x = 0.05 m in Fig. [5] Five time steps of 
t = 0.0001, 0.0005 , 0.001, 0.0015 and 0.002 s are presented. The potential along the crack is almost 
independent of the mesh-size, which is achieved by the special definition of a c in Eq. (fT2")) . Perpendicular 
to the crack, the potential has its maximum at the crack and decreases away from the crack (Fig. [HJ. 
This imbibition process is described independently of the mesh size. 

3.3 Fracture and flow in a random mesh 

The last example involves crack propagation and flow in a random mesh. In the previous example 
the elements in the middle of the specimen were aligned along the crack path, which was useful to 
evaluate the potential along and perpendicular to the crack. However, lattices for the analysis of fracture 
processes are usually random , since regular arrangements of lattice elements influence the direction of 
crack propagation [18| . In the present example, the mechanical lattice model is used to analyse splitting 
cracks due to expansion of an inclusion for three random meshes with varying element sizes. The specimen 
geometry and loading setup was chosen according to corrosion experiments in [l| , i.e. the circular inclusion 
corresponds to a steel reinforcing bar cross-section. The coarse mesh for structural analysis is shown in 
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Figure 5: Meshes and loading set-up for (a) structural and (b) transport analysis for the medium density 
lattice with elements aligned along the potential crack path. 




Figure 6: Results of the structural analysis: Load-displacement curves for three element sizes. The open 
circle marks the stage for which the flow analyses are performed in the second stage of the analysis. 
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Figure 7: Potential in x-direction along the crack at y = 0.05 m for several time steps. 




Figure 8: Potential in y-direction perpendicular to the crack at x = 0.05 m for several time steps. 
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Figure 9: Random lattices with inclusions for (a) structural and (b) transport analysis for d m = 0.002 m 
in the refined regions of the specimens. 



Fig. Within the square specimen of length L = 0.15 m, an inclusion of diameter d a = 16 mm is 
located with its center at x = 0.075 m and y = 0.12 m. The expansion of the inclusion is described by 
eigen-displacement ut subjected to elements, which cross the boundary of the inclusion. In these elements 
the normal component of the eigenstrain in Eq. ^ is £t = Ut/H. Thus, the eigenstrain depends on the 
lattice clement size, whereas the displacement is independent. A displacement of ut = 0.0665 mm at 
these elements is applied incrementally. The eigenstrain is applied uniformly along the circumference of 
the reinforcement bar. However, the total strain of the lattice elements crossing the circumference of the 
reinforcement bar is nonuniform, since the mechanical part of the strain depends on the stiffness of the 
surrounding material. The crack opening was evaluated as the relative displacement in x-direction of 
points A and B at the top of the specimen over a length l c = 0.016 m (Fig. [S^,). The model parameters 
were the same as for the previous example. In the refined area, the three lattices were generated with 
d m = 0.002, 0.001 and 0.0005 m. The results are presented in the form of inclusion expansion ut versus 
the crack opening in Figure 1101 The initial part of the expansion-crack opening curve is nearly mesh- 
independent. In a later stage, the results for the three meshes differ. However, the difference does not 
appear to be due the element size, since the medium mesh overestimates the crack opening obtained from 
the fine and coarse mesh. 

In the second step of this example, a flow analysis was performed. The nodes at the top, left and right 
boundaries were subjected to a constant potential of 9 = 1. All other nodes within the specimen had 
an initial value of 9 — 0. The diffusivity a was determined according to Eq. (jTTJ) . The results of the 
potential along the x-direction for y = 0.135 m is shown in Figure [TT] for t = 0.00001 s. Similarly, as in 
the second example, the potential has its maximum at the location of the crack, which is not positioned 
at the same cc-coordinate since the random lattices are used. Nevertheless, the shape and magnitude of 
the potentials along this section of the specimen are independent of the element size. 



4 Conclusions 



In the present work, a lattice approach to model flow in cracked media is presented. The approach couples 
mechanical loading to flow analysis by relating the diffusivity of conduit elements to the equivalent crack 
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Figure 10: Results of the structural analysis: Crack opening versus eigendisplacement ut for three mesh 
sizes. The open circles mark stages at which the flow analysis is performed. 
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Figure 11: Results of the flow analysis: Potential 6 along the ^-direction (y = 0.135 m) in the refined 
region of the mesh for three element sizes at a eigendisplacement of ut = 0.0665 mm and t = 0.00001 s. 
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opening of mechanical lattice elements. Voronoi and Delaunay tessellations are used to define clement 
connectivities for the flow and mechanical lattices, respectively. The work resulted in the following 
conclusions: 

• The lattice of conduit elements with cross-sections obtained from the Delaunay triangulation results 
in an accurate description of stationary flow fields for uncracked homogenous materials. 

• The proposed coupling of mechanical loading with flow analysis results in a mesh-independent 
description of load-displacement curves and flow fields for lattice elements aligned aligned along 
cracks. 

• For random lattices, the position of cracks depends on the arrangement of lattice elements. However, 
the crack openings obtained are independent of the lattice size. Furthermore, flow fields for the 
cracked material can be described mesh-independently. 

The present lattice model is capable of describing the coupling of fracture and flow in concrete independent 
of the element size. The interaction of discrete and continuous flow, which can be problematic in hybrid 
discrete continuum approaches, is circumvented by modelling both flows by means of the same set of 
lattice elements, in which the discrete flow is smeared out over the width of the elements. Future work 
will concern the extension of this lattice approach to 3D and its application to the modelling of corrosion 
induced deterioration. 
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